Stability of the shell structure in 2D quantum dots 
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We study the effects of external impurities on the shell structure in semiconductor quantum 
dots by using a fast response-function method for solving the Kohn-Sham equations. We perform 
statistics of the addition energies up to 20 interacting electrons. The results show that the shell 
structure is generally preserved even if effects of high disorder are clear. The Coulomb interaction 
and the variation in ground-state spins have a strong effect on the addition-energy distributions, 
which in the noninteracting single-electron picture correspond to level statistics showing mixtures 
of Poisson and Wigner forms. 



I. INTRODUCTION 

The functionality of nanoelectronic devices strongly 
depends on the quality of the samples produced in the 
fabrication process. In semiconductor quantum dots 
(QD)/ there may exist external impurities, which have 
considerable effects on the measurable quantities such as 
the addition energy and the chemical potential. Gener- 
ally, the physical properties of QD's can be determined 
by the shape of the confining potential restricting the 
electrons in the semiconducting material. 

In a vertical QD formed in a double-barrier het- 
erostructure, the measured addition energies (cf. elec- 
tron affinities) correspond well to the shell structure of a 
two-dimensional parabolic (harmonic) model potential^ 
Remarkably good agreement between the experiment and 
theory has been obtained by using this model and taking 
the electron-electron interactions into account by employ- 
ing, e.g., the spin-density-functional theory (SDFT)ii4 
There are, however, always sample-specific variations^ 
such that in QD's with more than only a few electrons, 
the agreement between the experiments and theory be- 
comes generally worse, presumably due to non-symmetric 
deviations of the actual confining potential. 

Deviations from the circular symmetry can be induced 
by impurity particles migrated through the hcterostruc- 
ture layers. States bound to hydrogenic impurities were 
found in the single-electron tunneling experiment by 
Ashoori et al& Recently, a clear deformation in the trans- 
port spectrum of a resonant tunneling device was shown 
to be due to an ionized single or double-acceptor near 
the QDjZ" For the time being, there is no experimental 
data showing systematic impurity effects, but the even- 
tual serial production of QD structures can lead to a 
considerable fraction of distorted samples. 

Theoretical research of QD's containing impurities has 
been mostly based on either single-impurity systems 
studied with exact diagonalization'* and the SDFT^ or 
on QD's with a relatively high impurity densityiSiiS In the 
latter case, Hirose et al.^ have performed SDFT studies 
on the energy fluctuations and the effects of interactions, 
as well as on the role of spin>ii These systems have ex- 
perimental realization in disordered (no shell structure) 



QD's, where the fluctuations in the Coulomb blockade 
peak spacings show interesting statistics beyond the con- 
stant interaction model combined with the random ma- 
trix theorypi^iiSi 

In this paper we focus on the stability of the shell struc- 
ture in primarily parabolic few-electron QD's. Particu- 
larly, the effects of external impurities on the addition en- 
ergies are studied up to the high-disorder limit. The two 
main extensions to previous studies are (i) the use of an 
impurity model based on recent experimental findings^ 
and (ii) the computation of the complete addition energy 
spectrum up to 20 electrons with both a single impurity 
at an arbitrary position, and finally with an ensemble of 
one thousand impurity configurations. The calculations 
are performed using the SDFT with a remarkably fast 
response-function iteration method for the Kohn-Sham 
(KS) equationSfi^ The resulting spectra demonstrate that 
the shell structure is generally preserved, although in the 
mid-shell regime there can be drastic changes in the ener- 
getics. The statistics of the noninteracting level spacings 
can be characterized by a mixture of Poisson and Wigner 
distributions. Those shapes are strongly affected by the 
Coulomb interaction and by the ground-state spins. 

The organization of this paper is as follows. First in 
Sec. ^ we define the model system and the parameters 
for the external impurities. In Sec. IIIII we introduce the 
response-function scheme for the density-functional cal- 
culations. The results for a single repulsive impurity are 
presented in Sec. IIVI Finally the addition-energy distri- 
butions from one thousand configurations are shown in 
Sec. |V] The paper is summarized in Sec. IVII 



II. MODEL 

We model the quantum dot by using the conventional 
effective-mass approximation and choose the material pa- 
rameters for GaAs, i.e., the effective mass m* = 0.067 me 
and the dielectric constant e = 12.7. We apply a strictly 
two-dimensional model by assuming a negligible degree of 
freedom for electrons in the vertical direction. The many- 
electron Hamiltonian is given in the effective atomic 
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where defines the strength of the external confine- 
ment, which is assumed to be parabolic in shape. The 
Coulomb potential caused by the Ni^-^p impurity particles 
located in the vicinity of the QD and having charges \qk\, 
is written as 
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where Rfc, and dk define the lateral and vertical positions 
of the particles. This model has been shown to result in 
a good agreement with the measured transport spectrum 
of a vertical QD.^ 

In the first calculations, we include a single impurity 
with a negative unit charge on the QD plane {d = 0) and 
vary only the lateral position R. Further, when perform- 
ing a statistical analysis, we assume that there exist 1 — 10 
randomly distributed impurity particles migrated in the 
tunneling barrier of the QD. We choose the ranges for the 
impurity positions as < Rk < 10 and < dk < ^0,3- 
We set Qk = —1 which is in a qualitative accordance with 
the results from the tunneling experiments^ In addition, 
the minimum spacing between two impurities is expected 
to be 0.1 a*g which is of the order of the lattice constant 
for the layer material (e.g., AlGaAs). 



III. METHOD 

The numerical procedure for solving the electronic 
properties of the system defined above is based on 
the spin-density- functional theory (SDFT). In the 
self-consistent Kohn-Sham formulation, the effective 
Schrodinger equation to be solved iteratively is written 
as 



V2 

— + ^Ks(r) 



(3) 



where the Kohn-Sham potential Vks is a sum of the 
external potential defined above, the Hartree potential, 
and the exchange-correlation potential given as V^c(r) — 
SExc[p,S,]/Sp'^{r). Here p is the electron density, a de- 
notes the spin index, and ^(r) = [p^(r) — p^{r)]/p{r) is 
the spin polarization. For E^c we use the local spin- 
density approximation based on the functional provided 
by Attaccalite and co-workers 

A detailed description of the following algorithms can 
be found in Refs. IT^It^IiM We start with solving for 
the lowest n solutions of the eigenvalue problem Q by 
applying the evolution operator. 



repeatedly to a set of states {tpj, 1 < J < n}, and or- 
thogonalizing the states after every step. Instead of the 
commonly used second-order factorization in combina- 
tion with the Gram-Schmidt orthogonalization, we use 
the fourth-order factorization for the evolution opera- 
given by 
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For the orthonormalization, we diagonalize the matrix of 
overlap integrals and construct from these a new set of 
orthonormal states. Depending on the physical system 
this method can be faster by up to a factor of 100 in 
comparison to the second-order factorization. 

A widespread problem of the density-functional calcu- 
lations is the large number of iterations required to obtain 
a self-consistent solution. Namely, to maintain the sta- 
bility of the iteration process, it is necessary to keep the 
mixing parameter(s) small, which can lead to thousands 
of iterations. We overcome this problem with a procedure 
having its roots in the Hohenberg-Kohn theorem which 
states that the one-body density p(r) is the only truly in- 
dependent variable. Therefore an algorithm which solves 
for p(r) directly by applying a Newton-Raphson proce- 
dure is used. We define 
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as the density difference between two self-consistent it- 
erations. Here n"' is the occupation factor, <I>g are the 
orthogonalizcd solutions of Eq. Q , and p"" (r) is the den- 
sity used for the calculation of Vks- The sum is over all 
occupied (hole: h) states. Then the density correction 
Sp'^{r) is determined by a linear equation 
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Here d is the dimension of the system and e°''°' is the 
static dielectric function of a non-uniform electron gas^ 
It contains the zero-frequency Lindhard function and 
Vpl'({r,r') ^ 6V-£s(j)/Sp'^' (r'). To avoid the calcu- 
lation of unoccupied (particle: p) states we seek for 
an approximation for the static response function that 
only needs the calculation of occupied states. For that 
purpose, we recall that linear response theory can be 
derived^^ from an action principle for excitations of the 
form \i>{t)) = e:>ip{J2ph Cph{t)alah) \cj)o), where \cj)(i) is the 
ground state, and Cph{i) are particle-hole amplitudes. If 
we assume that the particle-hole amplitudes are matrix 
elements of a local one-body operator c'^(r), we end up 
with Feynman's theory of collective excitationsi^ In this 
approximation, we can rewrite Eq. Q as 



= 2Y^S^p^V;-_i^ ^p^\ (8) 
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FIG. 1: Addition energies for clean and impurity-containing 
QD's with the confinement strengths of luq = 5 meV (a) and 
3 meV (b). 



where now 




FIG. 2: Noninteracting single-electron energies for a quantum 
dot containing a single repulsive impurity particle on the dot 
plane (d — 0) and at R = . . .8a*g from the center. 



and 
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which is proportional to the static structure function of 
the noninteracting system. Above, the asterisk stands for 
the convolution integral. With these manipulations, we 
have rewritten the response-iteration equation in a form 
that requires only the calculation of the occupied states. 
Since the multiplication of the left-hand side of Eq. © 
requires only vector- vector operations, the equation can 
be solved with the conjugate-gradient method. 



IV. SINGLE IMPURITY 

Figurenshows the addition energies up to 20 electrons 
for two QD's with luq — ^ and 3 meV. A single impurity 
particle with g = — 1 is located on the QD plane {d = 0) 
at i? = 2lo, where lo = ^JlJmFuJo, giving R = 3.00 
and 3.87 a|j for ujq — 5 and 3 meV, respectively. The 
addition energies of the corresponding clean (no impu- 
rity) QD's are also presented. They qualitatively agree 
with the previous SDFT results by Reimann et alf^ and 
Hirose and Wingreen^ for cjo 3 meV [Fig. ^b)]. In- 
creasing the confinement strength to 5 meV leads to the 
disappearance of the even peaks at = 14 ... 18, but 
the ground-state spins remain the same, having the well- 
known sequence determined by Hund's rule^ 

In the presence of the impurity, the peaks at magic 
electron numbers {N = 2, 6, 12, 20) are preserved. In con- 
trast, the peaks at iV = 4 and 9 which in the clean case 



correspond to aligned spins in the 2nd and 3rd shells, are 
strongly diminished and the ground-state spins are S = 
and 1/2, respectively. Generally, the inclusion of a single 
impurity worsens the agreement between the calculated 
addition-energy spectrum and the experimental result by 
Tarucha and co-workers^ 

Next we study how the energetics depends on the lo- 
cation of the impurity particle on the QD plane. The re- 
sults presented below have been calculated using the fixed 
confinement strength of = 5 meV. Figure |21 shows the 
noninteracting single-electron energies as a function of 
R. At large i?, the energies approach the shells of a 
clean QD with e„; — (2n -I- |/| -I- l)woj whereas R — Q 
corresponds to the eigenenergies of a finite-width quan- 
tum ring^'^ with several doubly-degenerated shells. The 
mid-i? range is characterized by many avoided crossings 
between the states that a degenerate in the absence of 
the impurity. 

We performed the SDFT calculations for iV = 2 ... 21 
and R — Q . . . 7.5 a*^ with a step size of 0.1 a*^. For each 
N and i?, we calculated several states to find the correct 
ground-state energy and spin. The plot of the addition 
energies is shown in Fig. 13 The values indicated by dif- 
ferent colors in the figure vary between 0.22 Ha* (dark 
blue) and 0.65 Ha* (dark red). The dots with magic elec- 
tron numbers remain particularly stable and paramag- 
netic {S = 0) despite the impurity particle located at 
different radius from the center. The most apparent ex- 
ception is TV = 12 at small i?, showing a shift to iV = 10 
as a stable configuration. This corresponds to a state 
with filled shells and = in the quantum-ring-like ge- 
ometry, as seen in Fig. El QD's with TV = 4, 8, 14, and 18 
have partially occupied highest shells in the both limits, 
leading to S* = 1 due to Hund's rule. The constant S = 1 
ground state of = 8 and 14 can be understood from 
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FIG. 3: (color online) Addition energies (scale from dark blue 
to dark red) for a quantum dot containing a single repulsive 
impurity particle at 7? = . . . 7.5 ag from the dot center. The 
filled circles denote the changes in the ground-state spin. The 
solid, dashed, and dotted lines correspond to S = or 1/2, 
S = 1 or 3/2, and 5* = 2 or 5/2, respectively. 




FIG. 4: Electron densities of a 12-electron quantum dot with 
an impurity particle located at 7? = 0, 1, 3, and 6 a*g from the 
dot center, respectively. 



the small energy gap between the 4th and 5th, as well as 
between the 7th and 8th eigenstates throughout the R 
regime (Fig. 13). 

In Fig. 0] we show examples of the electron densities for 
a 12-electron QD with different impurity locations. The 
insertion of the impurity particle at the center {R = 0) 
leads to a typical quantum-ring-like solution with the 
density distribution on the edge of the QD. As the impu- 
rity is located off the center, the density deforms grad- 
ually through a slightly localized shape with six peaks 
to a clean configuration where the circular symmetry is 
retained. 

We also find many solutions, always corresponding 



to half-filled shells, where the spin symmetry is bro- 
ken. These so-called spin-density waves (SDW)^ have 
been criticized by Hirose and Wingreen^, arguing that 
the solutions are mixtures of different (total) spin states, 
and thus artifacts of the SDFT. Recently, the claim was 
shown true explicitly in a rectangular four-electron QD 
by employing the SDFT and the exact diagonalization.— 
Nevertheless, the total energy of the SDW solution is only 
slightly smaller than the exact energy and in any case it 
represents an improvement over the corresponding single- 
configurational DFT (spin-compensated) energy. Hence, 
the SDW solutions are included as such in the energetics 
considered in this paper. 



V. RANDOM IMPURITIES 

Then we include a random number (1-10) of negatively 
charged impurities randomly in the expected tunneling 
barrier below the QD (see Sec. |n] for the parameters). 
We computed 1000 impurity configurations such that for 
each configuration the number of electrons in the QD 
was TV = 1 ... 21. To determine the groud-state spins in 
the many-electron problem, we calculated all the relevant 
{S < 2) spin states for each N. 

A. Single-electron picture 

We begin with the level statistics of a QD with nonin- 
teracting electrons. Now the electrons are affected only 
by the confining and impurity potentials. This corre- 
sponds to the single-electron picture where the addition 
energy is simply equal to the level spacing, Ao(A^) = 
ejv/2-i-i ~ eAf/2j where the divisor of two follows from the 
spin degeneracy. Fig. [51 shows examples of the (nonin- 
teracting) addition-energy distributions as = 4 and 
12. As the average impurity density is relatively low, 
there are several configurations that are close to the clean 
case, leading to a Poisson-like level distribution. On the 
other hand, the strongly distorted configurations lead to 
Wigner distribution. Thus, the sum of these subgroups 
give the observed double-peak shape which is pronounced 
in the 12-electron case [Fig. I^Jb)]. The active region of 
the QD is then larger than that for = 4, leading to a 
higher number of distorted configurations since the im- 
purity density is kept constant. In both cases, we have 
divided the configurations into these two parts in such 
a way that QD's with a total energy higher than a cer- 
tain, qualitatively chosen critical value are considered as 
"distorted" and the remaining configurations as "clean" . 
The fits corresponding to the sums of the Poisson and 
Wigner functions and Aq is plotted in Fig. [S] The distri- 
butions are mirror images of each other, since in the clean 
QD A^ = 12 is a fiUed-shell case with Ao(12) - 0.45 Ha*, 
whereas Ao(4) = 0. 

For QD's, the level-spacing distributions are typically 
plotted by keeping the QD parameters as the external 
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FIG. 5: Noninteracting addition-energy distributions of 1000 
impurity configurations corresponding to A'' = 4 (a) and 
A'^ = 12 (b). The approximative contributions of the Poisson 
and Wigner-type distributions are also shown. The energy 
binwidth is 0.01 Ha*. 



confinement or the magnetic field constant and consid- 
ering high number of eigenstates within a certain energy 
range. The integrability of the system can then be con- 
nected to the level spacings such that the Wigner distri- 
bution is interpreted as a signature of quantum chaosi2& 
Our case is qualitatively similar in the sense that the 
varied parameter is the impurity configuration instead of 
the energy, and high distortion leads to the Wigner dis- 
tribution as shown by Hirose and co-workers.^ We sug- 
gest, however, that in isolated few-electron QD's the av- 
erage impurity density is considerably smaller than in 
their model, leading to a significant contribution of the 
Poisson-type level-spacing distribution as visualized in 
Fig. |31 We note here that Hirose et al.^ defined the non- 
interacting picture such that the electrons do not interact 
with the impurities, while we consider the eigenenergies 
of a single electron in the presence of the impurity po- 
tential. 



B. Interacting electrons 

We examine the corresponding many-electron problem 
using the same 1000 randomly determined impurity con- 
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FIG. 6: Addition-energy distributions of 1000 different impu- 
rity configurations in a quantum dot of 1-20 electrons. All the 
possibilities for the ground-state spins are taken into account. 



figurations. Taking the different possibilities for each 
ground-state spin into account, we performed ~ 70000 
self-consistent SDFT calculations. FigureElshows the ad- 
dition energy spectrum, where the distribution for each 
N is plotted in grayscale. The mean values of the distri- 
butions are shown as pluses, and the crosses correspond 
to the clean, parabolic QD [repeated from Fig. H^a)]. 
Particularly with small N, the result of the clean case 
matches rather accurately with the maxima of the dis- 
tributions due to the relatively high fraction of non- 
distorted configurations. In the regime 13 < A^ < 19, 
there are clear deviations between the clean QD and the 
maxima, as visualized in the inset of Fig. El in more de- 
tail. The reason is the high degeneracy of different states 
in this regime and the increased tendency to distortion 
due to the relatively large active-dot region. 

In Fig.[7|we have a detailed look on the addition-energy 
distributions of A^ = 6 and 20 QD's in comparison with 
the corresponding values for noninteracting QD's defined 
in Sec. IV Al We have chosen these magic electron num- 
bers to exclude the spin effects. Namely, in determining 
A(6) and A(20), practically all impurity configurations 
give S = and 1/2 for these A^ and A^ ± 1, respectively. 
As expected, the Coulomb interaction shifts the distri- 
butions to higher energies, particularly in the iV = 6 
QD where the interaction strength is relatively larger. 
In general, the interactions enhance the average addi- 
tion energy of filled-shell QD's by a factor of 1.17 — 1.50 
(see Table U). In the mid-shell regime where the level- 
spacings fluctuate from degeneracy, the average nonin- 
teracting addition energies (Aq) are considerably lower, 
i.e., (A) -3... 7 (Ao). 

As seen in Fig. [7\ and Table Q interactions de- 
crease the fluctuation o f the addition energies defined 
as SA — ^ ((A — (A))^). The result is in contrast with 
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FIG. 7: Addition-energy distributions for A'' = 6 (a) and TV = 
20 (b) quantum dots in noninteracting (wliite) and interacting 
(dark gray) cases. 



TABLE I: Average noninteracting addition energies (Aq), 
their rms fluctuations SAo, the corresponding interacting case 
with (A) and SA, and the variation coefficient 5 A/ (A) for 
even A'^. The fractions of the ground-state spins are also 
shown, excluding A'^ = 16 quantum dot with Ps=2 = 0.28. 
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SA 


SA/ (A) 


Ps^o 
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2 


0.411 


0.0472 


0.615 


0.0521 


0.085 


1 





4 


0.070 


0.0595 


0.417 


0.0324 


0.078 


0.29 


0.71 


6 


0.353 


0.0777 


0.479 


0.0511 


0.107 


1 
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0.048 


0.0400 


0.327 


0.0193 


0.059 


0.11 


0.89 


10 


0.071 


0.0552 


0.294 


0.0191 


0.065 


0.37 


0.63 


12 


0.306 


0.0993 


0.383 


0.0428 


0.112 


1 





14 


0.042 


0.0339 


0.276 


0.0144 


0.052 


0.06 


0.94 


16 


0.045 


0.0348 


0.257 


0.0113 


0.044 


0.24 


0.48 


18 


0.070 


0.0517 


0.258 


0.0189 


0.073 


0.39 


0.61 


20 


0.271 


0.1119 


0.318 


0.0377 


0.118 


0.99 


0.01 



strongly disordered QD's^ and can be explained with the 
shell structure determined by the primarily symmetric 
(parabolic) confining potential: in the interacting system 
the impurities are screened by the electrons such that the 
QD is effectively more symmetric than in the correspond- 
ing noninteracting case, which leads to enhanced stabil- 
ity and smaller fluctuations. The variation coefficients 
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FIG. 8: Addition-energy distributions for two quantum dots 
of 4 and 14 interacting electrons, respectively. The results 
have beed calculated both by (i) taking all the possible 
ground-state spins into account (dark gray) and (ii) neglect- 
ing Hund's rule and thus setting S = and 1/2 for the ground 
states of A'^ = 4, 14 and A'^ = 3, 5, 13, 15, respectively (white). 



SA/ (A) show that the fluctuations are relatively largest 
for the magic electron numbers, since the conflgurations 
with high disorder have the strongest distorting effect on 
the corresponding flUed-shell structure. 

Figure |S1 shows the addition-energy distributions for 
A^ 4 and 14 QD's when the different ground-state spins 
arc included (dark grey) and neglected (white). With 
both N , there is (approximately) an equal possibility for 
S" = 1 (~ clean) and S = (distorted), leading to a 
larger width of the distribution. When TV = 4, the split- 
ting of the different spin states into two separate peaks is 
clear. Hence, the measured addition energy may depend 
strongly on the impurity conflguration of the particular 
sample. 

Fractions of 5 = and S — 1 ground states in the cal- 
culated conflgurations for even N are shown in Tabled 
Except a few results in TV = 12 QD and ~ 1 % of the 
configurations in A/ = 20 QD, the magic electron num- 
bers correspond to 5 = states. Otherwise, the S = 1 
states formed by Hund's rule are dominant. The TV = 16 
QD, which in the clean case has a half- filled shell with 
a four-fold degeneracy, has also a significant fraction of 
S — 2 states. 

We underline that due to the primarily parabolic con- 
fining potential and a relatively small number of impuri- 
ties, the shell structure is generally a dominating feature 
in the energetics. However, the features of disorder are 
clearly visible in the results: (i) A considerable fraction 
of the configurations shows Wigner shape in the level- 
spacing distribution (Fig. O; (ii) the relative fluctuation 
in (A) is largest for flUed-shell QD's; (iii) interactions 
bring Gaussian contribution to the distribution functions 
in the mid-shell regime (Fig.lSJ. The last feature is typ- 
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ical for chaotic and disordered QD's with considerable 
electron-electron interactions»i^ 



VI. SUMMARY 

To summarize, we have used a fast response-function 
algorithm within the spin-density-functional theory to 
study the stability of the shell structure in few-electron 
quantum dots. A single impurity near the quantum dot 
may cause dramatic changes in the addition energies, 
but the shell structure is generally preserved with high 
peaks at the magic electron numbers. The same effect 
was found using a model of one thousand random impu- 
rity configurations, even if the fluctuations are largest in 
the filled-shell regime. In the noninteracting picture, the 
addition-energy distributions are combinations of Pois- 



son and Wigner forms. The electron-electron interac- 
tions induce a shift to higher energies and considerably 
smaller fluctuations. In the mid-shell regime, changes in 
the ground-state spin lead to a strong broadening of the 
distribution. Hence, in a series of measurements there 
can be large fluctuations in those addition energies, and 
the result for a particular quantum dot strongly depends 
on the impurities in the sample. 
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